npj Systems Biology and Applications

High-dimensional population balance equations (PBEs) provide a natural framework for modeling heterogeneous cell populations, but their direct numerical solution becomes computationally prohibitive when the internal state space contains many molecular variables. We propose a hybrid mechanistic-machine learning framework for reducing and simulating PBEs defined over high-dimensional intracellular coordinates. The cell population is described by a number density n(x, t), where x [isin] [R]N represents gene and protein states associated with macrophage activation. A dynamics-preserving autoencoder maps this state space to a low-dimensional latent coordinate z [isin] [R]d, with d << N, while retaining key qualitative features of the underlying gene regulatory network, including attractor structure and multistability. Mechanistic information from the original regulatory dynamics is used to construct interpretable drift and diffusion terms for the reduced latent-space PBE. The reduced PBE is solved using a stochastic Lagrangian particle representation, in which particles evolve according to stochastic differential equations (SDEs) corresponding to the latent drift and diffusion fields. The resulting latent-space solution is subsequently decoded and propagated back into the original state space to recover physically interpretable cellular dynamics. We demonstrate the framework on macrophage polarization under cytokine-dependent regulation, including gene knockout perturbations. Overall, the proposed framework provides a computationally tractable and mechanistically interpretable route for integrating single-cell genomic data with population balance models of cell-state dynamics.

Optimizing dose and schedule remains a central challenge in oncology drug development, particularly for immunotherapies where fixed dosing regimens often fail to account for patient specific heterogeneity in tumor-immune dynamics. Here, we present a hybrid quantitative systems pharmacology-reinforcement learning-Monte Carlo Tree Search (QSP-RL-MCTS) framework for personalized immunotherapy dosing that formulates dose selection as a sequential decision-making problem. The approach integrates a mechanistic QSP model of prostate cancer immunotherapy, transcriptomics informed virtual patient populations and data driven AI system comprising reinforcement learning and Monte Carlo tree search. Reinforcement learning is used to learn adaptive generalized dosing policies that optimize treatment outcomes across the population, while Monte Carlo Tree Search provides forward-looking evaluation of RL predicted dosing trajectories to refine patient-specific decisions. On benchmarking against fixed dosing regimens of ipilimumab, the remission rate of the proposed model (95.2%) was comparable to the highest fixed dosing regimen of 10 mg/kg per dose while the median total dose (72 mg/kg) of the proposed model designed regimen was comparable to the lowest fixed dosing regimen of 3 mg/kg per dose. The model is generalizable across different dosing protocols and can be extended to predict optimal dose under different therapeutic scenarios. Analysis of the learned dosing trajectories enables stratification of patients into distinct response groups and identifies drug activity rate as the dominant determinant of long-term treatment outcome. These results demonstrate how mechanistically guided artificial intelligence can transform population-level dose optimization into patient-specific, biologically interpretable treatment strategies for precision immuno-oncology.

Analysis of circadian patterns in time-series data requires computational methods that can accommodate several factors, including variable sampling resolution, replicate number, and missing values. Most existing tools simplify rhythmicity to a strict dichotomy based solely on a single p-value threshold. This leads to a level of uncertainty that affects many biological targets. We developed CODAC (Circadian Oscillation Detection Analysis and Comparison), a framework that integrates nonlinear constrained optimization with a multicriteria rhythmicity classification scheme to evaluate rhythmic patterns without relying on a single statistical cutoff. This approach allows CODAC to identify and exclude medium-confidence rhythms rather than force them into a rhythmic/arrhythmic dichotomy. CODAC comprises four modules: (i) CODAC_single estimates rhythmicity within a single group; (ii) CODAC_flex extends this to identify distinct waveform types within one group; (iii) CODAC_compare performs pairwise comparisons across two or more groups to detect rhythmic or arrhythmic changes; and (iv) CODAC_multi handles more complex designs involving multiple-group comparisons. Using in silico simulations and public transcriptomic datasets, we show that CODAC performs comparably to established methods while providing additional flexibility for rhythm classification and comparison. Taken together, CODAC provides a flexible and open-source package for circadian timeseries analysis with automated visualization tools.

Systems biology aims to develop predictive models that connect molecular mechanisms to cellular behavior. Genome-scale metabolic models are among the most widely used frameworks for integrating stoichiometric, thermodynamic, and omics-derived information to predict feasible metabolic phenotypes. However, cellular metabolism operates on timescales governed by enzyme kinetics and by the relationship between metabolic fluxes and metabolite pool sizes. In steady-state metabolic models, this relationship can be expressed in terms of metabolite turnover rates, defined as flux-to-pool-size ratios that quantify how rapidly metabolite pools are renewed. As a result, physiologically consistent steady-state solutions should not only satisfy mass-balance and thermodynamic constraints but also exhibit turnover rates consistent with enzyme-mediated cellular dynamics. Current constraint-based approaches can admit many steady-state flux-concentration states that do not account for turnover rates, resulting in phenotypes incompatible with realistic metabolic dynamics, even when multiple types of data are imposed. Here, we present METEOR-K, an optimization framework that links steady-state metabolic fluxes to metabolite concentrations via turnover rate constraints to identify dynamically plausible flux-concentration reference states. Because these constraints reshape the feasible solution space, we also introduce turnover-rate-aware sampling strategies to efficiently explore the resulting feasible region. We applied METEOR-K to models of increasing scope and scale, including a reduced glycolysis pathway, anaerobic E. coli, and near-genome-scale ovarian cancer models. METEOR-K narrowed the admissible steady-state solution space, reduced uncertainty in feasible flux-concentration states, and improved local dynamic behavior. In nonlinear ODE simulations of bioreactor cultivation and drug-response scenarios, METEOR-K-derived states produced intracellular response times compatible with growth-supporting metabolic operation and perturbation recovery. Overall, these results establish metabolite turnover rates as scalable biophysical constraints that improve the physiological consistency of steady-state metabolic modeling. Because turnover rates encode flux-to-pool-size timescale constraints, METEOR-K moves part of physiological-consistency assessment upstream of kinetic parameterization, yielding better-suited flux-concentration reference states for kinetic modeling and dynamic prediction.

Breast-conserving surgery (BCS; lumpectomy) is widely used to treat early-stage breast cancer, yet the complexity and patient-specific variability of postoperative cavity remodeling make healing trajectories and physical outcomes difficult to predict. Although inflammatory and vascular processes are central to these outcomes, mathematical models of tissue healing do not capture their coupled interactions or calibrate them against experimental data. Here, we extend a computational model of breast cavity healing to incorporate coupled inflammatory and vascular dynamics, including angiogenesis, oxygen transport, and inflammatory cell activity. The model is calibrated using preclinical porcine lumpectomy histology and literature data. Model parameters are inferred using a multi-task Gaussian process surrogate within a Bayesian inference frame-work to align predictions with experimental observations and quantify uncertainty. Thus, this work provides a mechanistically grounded framework for inflammatory and vascular remodeling during breast cavity healing and provides a foundation for patient-specific prediction of healing and physical outcomes following lumpectomy.

Modern medicine implicitly assumes that physiological responses to intervention are predictably determined by administered treatments. However, physiological systems containing intrinsic delays between the detection of a stimulus and the biological response may violate this assumption. We investigate the human glucose-insulin system as described by the Ultradian model and mathematically demonstrate that clinically relevant forcing protocols-such as pulsatile insulin delivery and step-wise glucose infusion, both commonly used in intensive care units (ICUs)-can induce sustained temporal chaos that may hamper accurate prediction of the physiological response. If not accounted for, these chaotic dynamics could create difficulties in achieving optimal dosing and timing when administering glucose and insulin in clinical or home care settings. This phenomenon, termed delay-induced uncertainty (DIU), arises from the interaction between physiological delay, intrinsic shear near a limit cycle, and external forcing. Using the Ultradian glucose-insulin model, we compute top Lyapunov exponents to quantify predictability. Across a range of pulsatile and step-wise forcing regimes, including stochastic amplitudes drawn from Markov processes, we observe positive Lyapunov exponents, indicating sustained chaos. Our results suggest that delayed endocrine regulation may fundamentally limit the predictive value of the models used to develop glycemic management strategies, with implications for clinical protocols in the ICU.

Cells respond to their environment through protein networks often dysregulated in cancer, making dynamical modelling crucial. Limitations in experimental data and computational resources motivate coarse-graining methods to build low-dimensional descriptions. Yet classical approaches to coarse-grained modelling rely on strong assumptions, leaving it unclear when partial experimental observations support reduced descriptions of system dynamics. Here we show that symbolic regression (SR) provides a data-driven way to test whether, and how compactly, the dynamics of a signalling system coarse-grain over the measured variables, and, when they do, infers mechanistically interpretable models. In synthetic enzyme systems, SR recovers Michaelis-Menten kinetics for the two-step mechanism and under three-step extensions. As data quality is degraded, SR simplifies toward effective kinetic laws while preserving correct theoretical limits. Applied to published time-resolved ERK phosphorylation data, SR identifies compact phospho-ERK rate laws in selected cancer-relevant gene overexpression contexts, yielding interpretable kinetic effects. A sparse neural ODE baseline requires few inputs where SR succeeds, but on average more where it fails, indicating that, where a reduced model is learnable at all, SR failure is associated with more complex dynamics that a simple mathematical model cannot describe. Together, these findings establish symbolic regression as a way to test when a compact coarse-grained description is warranted, generating hypotheses where one holds and motivating potential new measurements where it does not.

More than 50 years ago, Haralick and co-authors proposed a family of gray-level co-occurrence statistics that became known as textural features. These features are widely used in image analysis, but their application to biological networks has remained limited because cellular networks are sparse, irregular graphs rather than regular pixel grids. This work presents a network-adapted version of Haralick texture analysis for generating pathway-level features from gene-level omics profiles. The resulting profiles reduce dimensionality and can be used as candidate predictors of anti-cancer drug response. Performance of these features is compared with original gene expression variables and with pathway features from network enrichment analysis (NEA), whose robustness has been demonstrated previously. Although technically simpler than NEA, Haralick features showed comparable sensitivity. More importantly, selected Haralick features were preserved between in vitro drug screens and clinical treatment-associated survival analyses, supporting their potential use for prioritizing robust pathway-level drug-response correlates.

Cancer cells rewire metabolism to support proliferation. Intriguingly, divergent metabolic choices are made to attain this common goal. Identifying the unique metabolic requirements for a specific cell has profound implications for cancer biology and precision medicine. Genome-scale metabolic models (GEMs) have emerged as powerful tools to systematically characterize, understand, and predict metabolism of cells and tissues. Despite being comprehensive, the current GEMs remain limited in their predictive power. Here, we present a new GEM of human cells, in silico Human Metabolic Essentiality (iHME), that significantly improves the prediction of metabolic dependencies at a reduced computational cost. Wse rationally downsized, curated, and corrected previous models to remove unsupported metabolic redundancies, which led to a slim model containing 4,377 reactions, 3,241 metabolites, and 1,825 genes. When used to reconstruct metabolic networks of 1,103 cancer cell lines, iHME recalled on average 84.6% of experimental essential genes, which is two-fold increase over previous models. Cholesterol biosynthesis was revealed to be the most reliably predicted pathway with alternative dependencies. Finally, we applied the model to reconstruct individualized networks and predict essential gene profiles for 8,384 patient tumor samples. Glucose transporter SLC2A1 (GLUT1) was identified as a context-specific dependency for head and neck cancers and ovarian cancer. Likewise, CDP-diacylglycerol synthase CDS2 was identified for skin cancer. Overall, iHME is a new genome-scale model for prediction of metabolic dependency at higher accuracy and computational efficiency.

Cancer cells rapidly rewire their metabolism, from efficient energy production toward anabolic processes, to sustain uncontrolled growth. Decoding such metabolic shifts is essential for uncovering novel therapeutic targets. To map systems-level metabolic changes across cancer types, we built context-specific genome-scale metabolic models for eight tissues (lung, thyroid, stomach, prostate, liver, kidney, colon, and breast) using gene expression data from The Cancer Genome Atlas (TCGA). Applying constraint-based modelling, we then identified differentially regulated pathways through flux enrichment analysis, revealing tissue-specific rewiring: branched chain amino acid metabolism was suppressed in breast cancer; sphingolipid metabolism was downregulated in colon, kidney, and thyroid but upregulated in breast. We further propose a model-driven pipeline to identify and characterise metabolic vulnerabilities. We first identify synthetic lethal reactions in normal tissues and their corresponding single lethal counterparts in cancers, thereby enabling the identification of metabolic "collateral lethal" reaction pairs for each cancer. Model-predicted collateral lethal gene pairs, including CMPK1-AK in colon, ALDOA-PGD in prostate, and SLC25A26-UQCRB in liver models, were supported through computational validation using DepMap data on gene essentiality. Subsequently, we show how to interpret metabolic rewiring in cancer tissues while accounting for any collateral lethal pairs. In summary, our results establish a systemic framework for decoding metabolic rewiring and synthetic lethal vulnerabilities in cancer.

Single-cell genome-scale metabolic models (scGEMs) enable characterization of metabolic heterogeneity underlying cellular states and phenotypes. However, methodological choices during scGEM construction can substantially alter model structure and predictions, challenging the reliability and comparability of resulting analyses. Here, we established a systematic benchmark to assess three key construction factors: data preprocessing method, model extraction method (MEM) and gene expression threshold. We evaluated 26 strategies representing different combinations of these factors across nine scRNA-seq datasets in three dimensions: accuracy, sensitivity to expression perturbation and computational feasibility. We found that MEM had the greatest influence on most accuracy metrics, data preprocessing strongly influenced the discrimination of cellular identities, and expression threshold balanced model completeness and cellular specificity. These findings indicate that strategy performance varies across evaluation criteria. Our benchmark provides practical guidance for strategy selection and an empirical basis for standardized evaluation and future scGEM method development.

A global pandemic requires accelerated strategies to mitigate public healthcare risks and preserve societal stability, either by repurposing existing drugs or by rapidly developing novel drug candidates. Computational platforms have accelerated candidate identification for both routes. A critical gap is the translation from in vitro potency to predicted clinical efficacy, which determines whether a prioritized candidate can achieve therapeutic effect under realistic dosing. To address this, we developed the Pandemic Pharmacology Platform for COMPUTational Evaluation of anti-infectives (PPP-COMPUTE), a comprehensive pharmacokinetic/pharmacodynamic (PK/PD) simulation platform designed to support evaluation and prioritization of drug candidates and clinical trial design during a pandemic. PPP-COMPUTE integrates experimentally derived anti-infective potency data (e.g., EC50) with PK/PD modeling to evaluate whether clinical dosing regimens can achieve sufficient exposure for therapeutic efficacy in patients. The platform incorporates mechanism-based dynamic models with a focus on viral pathogens to simulate time-dependent viral load trajectories under various treatment scenarios, enabling quantitative assessment of antiviral response and optimization of dosing strategies. Probability of target attainment analyses further support evaluation of regimen feasibility against predefined pharmacological targets. The clinical trial module generates simulated virological endpoints to evaluate candidate clinical study designs. PPP-COMPUTE is an accessible and quantitative framework that links PK and preclinical anti-infective potency data with predicted clinical benefit, thereby supporting rapid drug evaluation during future pandemics.

Skin growth and remodeling underlies health, disease, and treatments such as tissue expansion (TE). The mechanotransduction pathways in dermal fibroblasts are increasingly well characterized, and tissue-level growth has been described phenomenologically, but coupling between cell-level signaling and tissue-level growth remains poorly understood. We develop a dermal fibroblast signaling network through extensive literature data curation, comprising 151 reactions among 96 nodes. The inputs are mechanical stretch and eight ligands (TGF{beta}, PDGF, FGF, IL1, IL6, TNF, AngII, ET1); outputs of interest span ECM-enzymes (proMMP1/2/9, MMP1/2/9), ECM proteins (CImRNA, collagen I, fibronectin), and fibroblast activity (SMA, proliferation). Implemented as a logic-based ODE system, the network reproduces 82% of the calibration dataset and agrees with independent validation data. Sensitivity analysis reveals a tension-dependent regulation of signaling: at baseline tension, outputs are governed by many boosters (nodes that positively influence downstream targets) and one dominant brake, LATS1/2, whereas at high tension control consolidates and new, tension-specific regulators such as integrin (ITGB1) emerge. Multiple pathway axes converge on a few central regulators, producing pronounced crosstalk, most notably between TGF{beta} and mechanical tension. Finally, linking the collagen outputs to a tissue-level growth formulation yields a bidirectional mechanical-biochemical coupling that reproduces tension-induced skin growth measured in a porcine TE model. This framework establishes a comprehensively calibrated dermal fibroblast signaling network coupled to tissue-level growth, opening opportunities for targeted TE interventions.

Neutrophil swarming is a critical immune response in mammals and fish, in which neutrophils are recruited to inflammatory sites where they coordinate into a swarm that neutralizes pathogens. While excessive swarming can drive prolonged inflammation, a quantitative understanding of swarming dynamics remains limited. We developed a one-dimensional radial reaction-diffusion model of neutrophil swarming with two kinetic parameters, in order to capture the self-limiting swarming dynamics in both murine and human neutrophils in response to different inflammatory stimulus sizes. To ensure that the inverse problem is well-posed, we first performed sensitivity and identifiability analyses. We then developed a physics-informed neural network (PINN) to infer the key parameters governing swarm expansion and self-limitation. To account for uncertainty in noisy experimental measurements, we further extended this framework to a Bayesian PINN (B-PINN), which provides credible intervals for the inferred parameters. Both models were validated against synthetic data generated by numerical simulation and subsequently applied to in vitro experimental data from human and murine neutrophils in response to three bioparticle cluster sizes. The PINN-inferred dynamics show that larger bioparticle clusters are associated with greater cumulative recruitment and larger swarms in both species. The models further reveal species-specific differences in both the amplitude of initial recruitment and the timescale on which it self-limits. Additionally, the B-PINN posterior distributions quantify uncertainty in these species- and cluster size-dependent trends and identify where additional measurements would be most informative. To our knowledge, this is the first application of physics-informed machine learning to model neutrophil swarming dynamics. This framework provides a starting point for systematically comparing recruitment dynamics between human and murine neutrophils and offers guidance for future experimental design.

Isotopic Metabolic Flux Analysis (I-MFA) is a standard approach for estimating intracellular metabolic fluxes. I-MFA infers fluxes by comparing simulated and measured metabolite isotopologue distributions (MIDs) of metabolites from isotope labeling experiments. MIDs represent fractional abundances that strictly sum to one for any given metabolite, thus they are inherently compositional data. However, state-of-the-art estimation approaches rely on calculating standard Euclidean distances between MIDs in a non-compositional paradigm, introducing a systemic bias. To resolve this, our study proposes compositional I-MFA. We demonstrate how to construct a meaningful orthonormal basis for MIDs via ordered sequential binary partitioning, which can be used to perform isometric log-ratio (ILR) transformation. As a minimal change to existing I-MFA workflows, we suggest estimating fluxes by minimizing Euclidean distances between ILR-transformed MIDs. We validated this framework against traditional methods using both a toy model and a biologically realistic model, evaluating point estimates, sensitivity across varied true fluxes, and confidence intervals. In the two examples, compositional I-MFA consistently outperformed traditional approaches, reducing mean squared error of flux point estimates by an average of 42.6% and substantially narrowing confidence intervals. We conclude that compositional data analysis significantly improves I-MFA and can be implemented as a simple drop-in replacement for current pipelines. Graphical Abstract O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=156 SRC="FIGDIR/small/742769v1_ufig1.gif" ALT="Figure 1"> View larger version (34K): org.highwire.dtl.DTLVardef@1a53aa4org.highwire.dtl.DTLVardef@ad225aorg.highwire.dtl.DTLVardef@aa430eorg.highwire.dtl.DTLVardef@1880ca_HPS_FORMAT_FIGEXP M_FIG C_FIG HighlightsO_LINew compositional data approach improves metabolic flux estimation. C_LIO_LIThis data transformation requires minimal changes to existing workflows. C_LIO_LIThe new method reduced MSE of flux estimates by 42.6% in two examples tested. C_LIO_LIThe confidence intervals of the estimated fluxes were substantially narrowed. C_LIO_LIEstimation accuracy remained robust across a wide range of metabolic fluxes. C_LI

Computational modeling provides a powerful framework for in silico exploration of anti-cancer therapeutic targets and tumor response mechanisms. Oncogenic signaling pathways play a central role in tumor behavior and represent promising targets for personalized combination therapies. However, these pathways are complex, and although logic-based models are well suited for representing signaling dynamics, they are often constrained by model-specific data requirements, limited scalability, and time-consuming manual curation. Here, we introduce Functional Integration of Contextualized Omics for Unraveling regulatory dynamicS (FICUS), a framework that integrates omics-driven network contextualization with dynamic Boolean and logic-ODE modeling. FICUS enables automated, data-driven protein network inference and patient stratification, allowing shared signaling mechanisms to be identified across patient subgroups while preserving patient-specific dynamic responses. We applied FICUS to the SU2C-MARK lung cancer cohort and the The Cancer Genome Atlas kidney cancer cohort, demonstrating its utility for post-hoc analyses and downstream interrogation of dynamic tumor models. Overall, our results highlight the flexibility of FICUS in capturing heterogeneous signaling mechanisms across patient subgroups, addressing a key challenge in precision oncology.

BackgroundAcute stress responses are often reversible, whereas sustained stress can produce coordinated disruption across endocrine, metabolic, inflammatory, antioxidant, and neuroplastic pathways. The Tiered Stress Biochemistry Model (TSBM) is a hypothesis-generating ten-state ordinary differential-equation framework linking a stylized cortisol signal to noradrenergic drive, vitamin C, a phenomenological aldosterone/renin-angiotensin drive, magnesium, normalized BDNF-related and Nrf2-related states, inflammation, and tryptophan-kynurenine metabolism. MethodsFive prespecified scenarios (normal, acute, chronic, depression-like, and low-cortisol PTSD-like) were simulated for 168 hours. Analyses included local stability, output-specific sensitivity screening, structural ablation, a 400-draw Latin-hypercube scan over independent {+/-} 20% parameter ranges, uncertainty distributions for threshold-crossing times, a success-conditioned parameter-trade-off screen, and a wider scan in which hypothesized BDNF-feedback gain magnitudes varied log-uniformly from 0.1 to 10 times nominal. Parameters were classified as literature-derived, literature-constrained/model-defined, calibrated, or hypothesized. ResultsUnder the specified forcing assumptions, sustained stress produced coordinated changes across several pathways. In the depression-like scenario, removing slow stress-load accumulation increased day-7 BDNF-related activity from 47.2% to 82.6%, reduced inflammation from 14.63 to 2.15 arbitrary units, and lowered KYN/TRP from 0.173 to 0.057. Removing BDNF-dependent gain increased BDNF only to 49.1% and delayed KYN/TRP crossing by 2.5 hours. The stricter multi-output conclusion was retained in 91% of uncertainty draws. Median crossing times retained the nominal sequence, but the complete four-event order occurred in only 43% of all draws. ConclusionsWithin this reduced model, a shared slow stress-load process coordinates the high-exposure state, whereas BDNF-dependent feedback acts mainly as an amplifier. The model generates an experimentally testable staging hypothesis: under sustained high-exposure forcing, magnesium changes may precede later BDNF-related and KYN/TRP changes. Longitudinal studies are required to determine whether this sequence occurs biologically, whether it is reversible, and whether it has clinical relevance. The simulations are not diagnostic or treatment recommendations.

In mammalian cells, under normal circumstances, the p53 protein exhibits oscillatory dynamics in response to DNA damage and maintains the cells in a cell-cycle-arrested state. Intriguingly, some cells can escape this cell-cycle-arrested state even after prolonged DNA damage, and often undergo mitotic catastrophe. In this context, the precise role of p53 dynamics and its complex interplay with cell-cycle regulation remain poorly understood. Herein, by constructing a comprehensive network model, we have identified crucial crosstalk regulations between the p53 protein and key cell-cycle regulators that enable some cells to escape cell-cycle arrest during prolonged DNA damage. The model further illustrates a probable cellular mechanism underlying mitotic catastrophe and predicts ways to induce it in a therapeutically relevant manner.

Formalin-fixed paraffin-embedded (FFPE) tissue is the dominant clinical pathology resource yet whether it faithfully preserves organ signalling biology and supports directional regulatory analysis remains unquantified. We generated a phosphoproteome map from eight healthy rat organs, separating preservation effects from biological variation. Using mass spectrometry, we quantified 54,710 phosphosites on 5,994 proteins across receptors, kinase cascades and nuclear regulators. Organ-specific phosphosite signatures matched known physiological and proliferative states. Paired antagonistic phosphosites converted into "activating-minus-inhibitory" indices that quantified net tissue-specific pathway activity, while a "kinase-by-organ activity" matrix resolved functional hierarchies. Joint analysis with an external fresh-frozen phosphoproteome dataset yielded 58,631 phosphosites total, recovering 86% of the 28,888 sites detected in the frozen dataset. Organ identity explained over 92% of the total variance after batch correction, versus under 0.5% for preservation method. Per-organ phosphosite intensities agreed closely between preservation modes except in brain. This establishes that archived pathology tissue supports biologically faithful phosphoproteome analysis at organ, pathway, and site resolution, providing a framework for retrospective signalling studies in clinical archives. Graphical Abstract O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=105 SRC="FIGDIR/small/741173v1_ufig1.gif" ALT="Figure 1"> View larger version (44K): org.highwire.dtl.DTLVardef@7ec1forg.highwire.dtl.DTLVardef@1f22bcorg.highwire.dtl.DTLVardef@217f5forg.highwire.dtl.DTLVardef@1316b61_HPS_FORMAT_FIGEXP M_FIG C_FIG

We present DynaVelo, a generative neural ordinary differential equation model that learns the joint dynamics of gene expression and transcription factor (TF) motif activities in evolving cell systems using single-cell multiome with joint gene expression and chromatin accesibility readout. DynaVelo leverages partial RNA velocity information together with single-cell TF motif accessibility data to improve the modeling of cell state dynamics and identification of TF drivers. We show that DynaVelo recovers the complex and bifurcating in vivo dynamics of wildtype murine germinal center (GC) B cells and reveals how these cell dynamics change under loss-of-function mutations in epigenetic regulators Arid1a and Ctcf. DynaVelo resolves how TF motif activities evolve along latent time trajectories using analysis of training cells or through generated trajectories from the model. In silico perturbation analysis further enables DynaVelo to infer dynamic and cell-state-specific gene regulatory networks (GRNs), recovering many known TF-to-gene edges in the wildtype GC GRN and predicting those that are disrupted in mutants. Finally, in silico gene and TF perturbations allow both the prediction of cell dynamics under loss-of-function genetic mutations and the identification of TF perturbations to rescue loss-of-function dynamic and immunological phenotypes. This analysis predicted that Ctcf knockout would rescue Arid1a loss-of-function phenotype in the GC reaction and nominated Bcl6 and Stat3 as additional TFs whose knockout would rescue Arid1a loss. We validated these predictions in vivo using double heterozygous mutant mice, confirming rescue of the Arid1a dark zone phenotype in all cases and quantitatively assessing model predictions using multiome in Arid1aHet;CtcfHet double heterozygous mice. DynaVelo therefore provides a powerful new deep learning framework for modeling and perturbing dynamic cell systems by harnessing single-cell multiome data sets.